Assignment #3

# SEE modeldata package for new datasets
library(tidyverse)         # for graphing and data cleaning
library(tidymodels)        # for modeling
library(stacks)            # for stacking models
library(naniar)            # for examining missing values (NAs)
library(lubridate)         # for date manipulation
library(moderndive)        # for King County housing data
library(DALEX)             # for model interpretation  
library(DALEXtra)          # for extension of DALEX
library(patchwork)         # for combining plots nicely
library(dbplyr) 
library(scales)                # for SQL query "cheating" - part of tidyverse but needs to be loaded separately
library(mdsr)              # for accessing some databases - goes with Modern Data Science with R textbook
library(kableExtra)
library(RMySQL)            # for accessing MySQL databases
library(RSQLite)           # for accessing SQLite databases

#mapping
library(maps)              # for built-in maps
library(sf)                # for making maps using geom_sf
library(ggthemes)          # Lisa added - I like theme_map() for maps :)

#tidytext
library(tidytext)          # for text analysis, the tidy way!
library(textdata)          
library(reshape2)
library(wordcloud)         # for wordcloud
library(stopwords)

theme_set(theme_minimal()) # Lisa's favorite theme

When you finish the assignment, remove the # from the options chunk at the top, so that messages and warnings aren’t printed. If you are getting errors in your code, add error = TRUE so that the file knits. I would recommend not removing the # until you are completely finished.

Put it on GitHub!

From now on, GitHub should be part of your routine when doing assignments. I recommend making it part of your process anytime you are working in R, but I’ll make you show it’s part of your process for assignments.

Task: When you are finished with the assignment, post a link below to the GitHub repo for the assignment. If you want to post it to your personal website, that’s ok (not required). Make sure the link goes to a spot in the repo where I can easily find this assignment. For example, if you have a website with a blog and post the assignment as a blog post, link to the post’s folder in the repo. As an example, I’ve linked to my GitHub stacking material here.

https://github.com/jakebulling/HW3.git

Local Interpretable Machine Learning

You are going to use the King County house data and the same random forest model to predict log_price that I used in the tutorial.

data("house_prices")

# Create log_price and drop price variable
house_prices <- house_prices %>% 
  mutate(log_price = log(price, base = 10)) %>% 
  # make all integers numeric ... fixes prediction problem
  mutate(across(where(is.integer), as.numeric)) %>% 
  select(-price)


set.seed(327) #for reproducibility

# Randomly assigns 75% of the data to training.
house_split <- initial_split(house_prices, 
                             prop = .75)
house_training <- training(house_split)
house_testing <- testing(house_split)

ranger_recipe <- 
  recipe(formula = log_price ~ ., 
         data = house_training) %>% 
  step_date(date, 
            features = "month") %>% 
  # Make these evaluative variables, not included in modeling
  update_role(all_of(c("id",
                       "date")),
              new_role = "evaluative")

#define model
ranger_spec <- 
  rand_forest(mtry = 6, 
              min_n = 10, 
              trees = 200) %>% 
  set_mode("regression") %>% 
  set_engine("ranger")

#create workflow
ranger_workflow <- 
  workflow() %>% 
  add_recipe(ranger_recipe) %>% 
  add_model(ranger_spec) 

#fit the model
set.seed(712) # for reproducibility - random sampling in random forest choosing number of variables
ranger_fit <- ranger_workflow %>% 
  fit(house_training)

Tasks:

1. Choose 3 new observations and do the following for each observation:
   

• Construct a break-down plot using the default ordering. Interpret the resulting graph. Which variables contribute most to each observation’s prediction?

rf_explain <- 
  explain_tidymodels(
    model = ranger_fit,
    data = house_training %>% select(-log_price), 
    y = house_training %>%  pull(log_price),
    label = "rf"
  )

## Preparation of a new explainer is initiated
##   -> model label       :  rf 
##   -> data              :  16210  rows  20  cols 
##   -> data              :  tibble converted into a data.frame 
##   -> target variable   :  16210  values 
##   -> predict function  :  yhat.workflow  will be used (  default  )
##   -> predicted values  :  No value for predict function target column. (  default  )
##   -> model_info        :  package tidymodels , ver. 0.1.2 , task regression (  default  ) 
##   -> predicted values  :  numerical, min =  5.081823 , mean =  5.66507 , max =  6.726297  
##   -> residual function :  difference between y and yhat (  default  )
##   -> residuals         :  numerical, min =  -0.3147814 , mean =  0.0004295329 , max =  0.2270978  
##   A new explainer has been created! 

obs4 <- house_training %>%
  slice(4)
obs250 <- house_training %>%
  slice(250)
obs1000 <- house_training %>%
  slice(1000)
plot(predict_parts(explainer = rf_explain,
              new_observation = obs4,
              type = "break_down"))

plot(predict_parts(explainer = rf_explain,
              new_observation = obs250,
              type = "break_down"))

plot(predict_parts(explainer = rf_explain,
              new_observation = obs1000,
              type = "break_down"))

The latitude variables are the ones that have the largest effect on the prediction and square foot living is near the largest effect as well.

• Construct a SHAP graph and interpret it. Does it tell a similar story to the break-down plot?

rf_shap <-predict_parts(explainer = rf_explain,
                        new_observation = obs4,
                        type = "shap",
                        B = 10)

rf_shap_250 <-predict_parts(explainer = rf_explain,
                        new_observation = obs250,
                        type = "shap",
                        B = 10)
rf_shap_1000 <-predict_parts(explainer = rf_explain,
                        new_observation = obs1000,
                        type = "shap",
                        B = 10)

plot(rf_shap) 

plot(rf_shap_250) 

plot(rf_shap_1000)

These results are consistent with the results obtained earlier that latitude, longitude, and the square foot living are the largest contributors.

• Construct a LIME graph (follow my code carefully). How close is each original prediction to the prediction from the local model? Interpret the result. You can also try using fewer or more variables in the local model than I used in the example.

set.seed(494)
model_type.dalex_explainer <- DALEXtra::model_type.dalex_explainer
predict_model.dalex_explainer <- DALEXtra::predict_model.dalex_explainer

lime_rf <- predict_surrogate(explainer = rf_explain,
                             new_observation = obs4 %>%
                               select(-log_price), 
                             n_features = 5,
                             n_permutations = 1000,
                             type = "lime")
plot(lime_rf) + 
  labs(x = "Variable")

2. Describe how you would use the interpretable machine learning tools we’ve learned (both local and global) in future machine learning projects? How does each of them help you?

I would use the interpretable machine learning tools to help me discover which variables I have that are most important to determining our predictions. This could help us further study the relationships that we observe from the interpretable machine learning tools.

SQL

You will use the airlines data from the SQL database that I used in the example in the tutorial. Be sure to include the chunk to connect to the database here. And, when you are finished, disconnect. You may need to reconnect throughout as it times out after a while.

Tasks:

1. Create a SQL chunk and an equivalent R code chunk that does the following: for each airport (with its name, not code), year, and month find the total number of departing flights, the average distance of the flight, and the proportion of flights that arrived more than 20 minutes late. In the R code chunk, write this out to a dataset. (HINT: 1. start small! 2. you may want to do the R part first and use it to “cheat” into the SQL code).

• With the dataset you wrote out, create a graph that helps illustrate the “worst” airports in terms of late arrivals. You have some freedom in how you define worst and you may want to consider some of the other variables you computed. Do some theming to make your graph look glamorous (those of you who weren’t in my intro data science class this year may want to watch Will Chase’s Glamour of Graphics talk for inspiration).

con_air <- dbConnect_scidb("airlines")

flights_df <- tbl(con_air, "flights") 


airports_df <- tbl(con_air, "airports")



flight_airport_df <- 
  flights_df %>%
  filter(year == 2017) %>%
  mutate(faa = origin) %>%
  group_by(origin, month) %>%
  summarize(num_depart = n(),
             avg_distance = mean(distance),
             avg_delay = mean(arr_delay > 20)) %>%
  inner_join(airports_df, by = c("origin" = "faa"))


flight_airport_df %>%
  show_query()

## <SQL>
## SELECT `origin`, `month`, `num_depart`, `avg_distance`, `avg_delay`, `name`, `lat`, `lon`, `alt`, `tz`, `dst`, `city`, `country`
## FROM (SELECT `origin`, `month`, COUNT(*) AS `num_depart`, AVG(`distance`) AS `avg_distance`, AVG(`arr_delay` > 20.0) AS `avg_delay`
## FROM (SELECT `year`, `month`, `day`, `dep_time`, `sched_dep_time`, `dep_delay`, `arr_time`, `sched_arr_time`, `arr_delay`, `carrier`, `tailnum`, `flight`, `origin`, `dest`, `air_time`, `distance`, `cancelled`, `diverted`, `hour`, `minute`, `time_hour`, `origin` AS `faa`
## FROM `flights`
## WHERE (`year` = 2017.0)) `q01`
## GROUP BY `origin`, `month`) `LHS`
## INNER JOIN `airports` AS `RHS`
## ON (`LHS`.`origin` = `RHS`.`faa`)

flight_airport_rdf <- flight_airport_df %>%
  collect()

SELECT `origin`, `month`, `num_depart`, `avg_distance`, `avg_delay`, `name`, `lat`, `lon`, `alt`, `tz`, `dst`, `city`, `country`
FROM (SELECT `origin`, `month`, COUNT(*) AS `num_depart`, AVG(`distance`) AS `avg_distance`, AVG(`arr_delay` > 20.0) AS `avg_delay`
FROM (SELECT `year`, `month`, `day`, `dep_time`, `sched_dep_time`, `dep_delay`, `arr_time`, `sched_arr_time`, `arr_delay`, `carrier`, `tailnum`, `flight`, `origin`, `dest`, `air_time`, `distance`, `cancelled`, `diverted`, `hour`, `minute`, `time_hour`, `origin` AS `faa`
FROM `flights`
WHERE (`year` = 2017.0)) `q01`
GROUP BY `origin`, `month`) `LHS`
INNER JOIN `airports` AS `RHS`
ON (`LHS`.`origin` = `RHS`.`faa`)

Displaying records 1 - 10

origin	month	num_depart	avg_distance	avg_delay	name	lat	lon	alt	tz	dst	city	country
ABE	3	167	560.6766	0.1617	Lehigh Valley Intl	40.65208	-75.44081	393	-5	A	Allentown	United States
ABE	4	139	586.3525	0.1367	Lehigh Valley Intl	40.65208	-75.44081	393	-5	A	Allentown	United States
ABE	5	160	574.4938	0.1125	Lehigh Valley Intl	40.65208	-75.44081	393	-5	A	Allentown	United States
ABE	6	134	572.2313	0.1119	Lehigh Valley Intl	40.65208	-75.44081	393	-5	A	Allentown	United States
ABE	7	192	573.5104	0.1354	Lehigh Valley Intl	40.65208	-75.44081	393	-5	A	Allentown	United States
ABE	8	180	562.2167	0.1000	Lehigh Valley Intl	40.65208	-75.44081	393	-5	A	Allentown	United States
ABE	9	206	606.0874	0.0777	Lehigh Valley Intl	40.65208	-75.44081	393	-5	A	Allentown	United States
ABE	10	253	593.4704	0.1344	Lehigh Valley Intl	40.65208	-75.44081	393	-5	A	Allentown	United States
ABE	11	220	585.3227	0.0909	Lehigh Valley Intl	40.65208	-75.44081	393	-5	A	Allentown	United States
ABE	1	191	579.8796	0.1466	Lehigh Valley Intl	40.65208	-75.44081	393	-5	A	Allentown	United States

flight_airport_rdf %>%
  filter(num_depart >1000) %>%
  mutate(name = fct_reorder(name, avg_delay)) %>%
  ggplot() + 
  aes(x = avg_delay, 
      y = fct_reorder(name, avg_delay)) + 
  geom_col(fill = "lightblue")+ 
  scale_x_continuous(expand = c(0,0), 
                     labels = percent) + 
  labs(x = NULL,
       y = NULL,
       title = "Which airlines had the largest % of flights that \nwere more than 20 minutes late from 2010-2017?") +
  theme_minimal()

• Although your graph was truly inspirational, you’ve been requested to “boil it down to a few numbers.” Some people just don’t appreciate all that effort you put in. And, you need to use the already summarized data that you already pulled in from SQL. Create a table with 6 or fewer rows and 3 or fewer columns that summarizes which airport is the “worst” in terms of late arrivals. Be careful with your calculations. You may consider using the kable, kableExtra, or gt packages to make your table look truly spectacular.

table_df<-flight_airport_rdf %>%
  select(num_depart, avg_delay, name) %>%
  group_by(name)%>%
  summarise(prop_late_over_20 = mean(avg_delay), 
            num_depart = sum(num_depart)) %>%
  filter(num_depart > 1000) %>%
  arrange(desc(prop_late_over_20)) %>%
  head(6)

table_df %>%
  kbl(caption = "Worst Airports by Proportion Flights late over 20 minutes.(min 1000 flights)") %>%
  kable_classic(full_width = F, html_font = "Cambria")

Worst Airports by Proportion Flights late over 20 minutes.(min 1000 flights)

name	prop_late_over_20	num_depart
Key Field	0.2485000	1052
Newark Liberty Intl	0.2142000	115968
Arcata	0.2122833	1202
Charles M Schulz Sonoma Co	0.2084091	1111
Burlington Intl	0.2076250	3274
San Francisco Intl	0.2063500	174631

2. Come up with your own interesting question that data in the airlines database can help you answer. Write a SQL query and equivalent R code chunk to extract the data you need and create an elegant graph to help answer the question. Be sure to write down the question so it is clear.

Which airport was the the most popular destination?

most_delay_time <- 
  flights_df %>%
  filter(year == 2015) %>%
  group_by(dest) %>%
  summarize(flights_received = n()) %>%
  inner_join(airports_df, by = c("dest" = "faa")) %>%
  ungroup() %>%
  group_by(name)

most_delay_time_df <- most_delay_time %>%
  collect()

most_delay_time_df%>%
  arrange(desc(flights_received))

most_delay_time_df %>%
  filter(flights_received > 10000)%>%
  ggplot() + 
  aes(x = flights_received, y = fct_reorder(name, flights_received)) + 
  geom_col(fill = "navy blue") + 
  labs(x = NULL, 
       y = NULL, 
       title = "Which airport recieved the largest amount of flights in 2015? ") + 
  scale_x_continuous(labels = scales::comma) + 
  theme_minimal()

Function Friday

If you need to revisit the material, it is posted on the moodle page. I’ve tried to add all the necessary libraries to the top, but I may have missed something.

geom_sf() tasks:

Using the example from class that we presented as a baseline (or your own if you really want to be ambitious), try to add the following components to the map of the contiguous United States:

1. Change the color scheme of the map from the default blue (one option could be viridis).

states <- st_as_sf(maps::map("state", plot = FALSE, fill = TRUE))
states <- states %>%
  mutate(area = as.numeric(st_area(states)))


ggplot(data = states) +
  geom_sf(aes(fill = area), color = "black") +
  coord_sf(xlim = c(-127, -63), ylim = c(24, 51), expand = FALSE) +
  scale_fill_gradient(low = "white", high = "purple") +
  theme_map()

2. Add a dot (or any symbol you want) to the centroid of each state.

?stat_sf_coordinates
ggplot(data = states) +
  geom_sf(aes(fill = area) , color = "black") +
  stat_sf_coordinates() +
  coord_sf(xlim = c(-127, -63), ylim = c(24, 51), expand = FALSE) +
  scale_fill_gradient(low = "white", high = "purple") +
  theme_map()

3. Add a layer onto the map with the counties.

counties <- st_as_sf(maps::map("county", plot = FALSE, fill = TRUE))
counties <- counties %>%
  mutate(area = as.numeric(st_area(counties)))



ggplot(data = states) +
    geom_sf(data = counties, aes(fill = area), color = "black") + 
  coord_sf() + 
  scale_fill_gradient(low = "light blue", high = "purple") + 
  theme_map()

4. Change the coordinates of the map to zoom in on your favorite state.

ggplot(data = states) +
    geom_sf(aes(fill = area), color = "black")+ 
  coord_sf(xlim = c(-109, -103), ylim = c(30, 37)) + 
  scale_fill_gradient(low = "light blue", high = "purple")  

Hint: https://www.r-spatial.org/r/2018/10/25/ggplot2-sf-2.html is a useful reference for some of the questions

tidytext tasks:

Now you will try using tidytext on a new dataset about Russian Troll tweets.

Read about the data

These are tweets from Twitter handles that are connected to the Internet Research Agency (IRA), a Russian “troll factory.” The majority of these tweets were posted from 2015-2017, but the datasets encompass tweets from February 2012 to May 2018.

Three of the main categories of troll tweet that we will be focusing on are Left Trolls, Right Trolls, and News Feed. Left Trolls usually pretend to be BLM activists, aiming to divide the democratic party (in this context, being pro-Bernie so that votes are taken away from Hillary). Right trolls imitate Trump supporters, and News Feed handles are “local news aggregators,” typically linking to legitimate news.

For our upcoming analyses, some important variables are:

• author (handle sending the tweet)
• content (text of the tweet)
• language (language of the tweet)
• publish_date (date and time the tweet was sent)

Variable documentation can be found on Github and a more detailed description of the dataset can be found in this fivethirtyeight article.

Because there are 12 datasets containing 2,973,371 tweets sent by 2,848 Twitter handles in total, we will be using three of these datasets (one from a Right troll, one from a Left troll, and one from a News Feed account).

1. Read in Troll Tweets Dataset - this takes a while. You can cache it so you don’t need to read it in again each time you knit. Be sure to remove the eval=FALSE!!!!

troll_tweets <- read_csv("https://raw.githubusercontent.com/fivethirtyeight/russian-troll-tweets/master/IRAhandle_tweets_12.csv")

2. Basic Data Cleaning and Exploration

1. Remove rows where the tweet was in a language other than English
2. Report the dimensions of the dataset
3. Create two or three basic exploratory plots of the data (ex. plot of the different locations from which tweets were posted, plot of the account category of a tweet)

#A
new_trolls <- troll_tweets %>%
  filter(language == "English")

#B

dim(new_trolls)

## [1] 175966     21

#C
ggplot(new_trolls) + 
  aes(x = account_category) + 
  geom_bar()

ggplot(new_trolls) + 
  aes(x = followers, y = following) + 
  geom_point()

3. Unnest Tokens

We want each row to represent a word from a tweet, rather than an entire tweet. Be sure to remove the eval=FALSE!!!!

troll_tweets_untoken <- new_trolls %>%
  unnest_tokens(word, content)

troll_tweets_untoken

4. Remove stopwords. Be sure to remove the eval=FALSE!!!!

#get rid of stopwords (the, and, etc.)
troll_tweets_cleaned <- troll_tweets_untoken %>%
  anti_join(stop_words)

Take a look at the troll_tweets_cleaned dataset. Are there any other words/letters/numbers that we want to eliminate that weren’t taken care of by stop_words? Be sure to remove the eval=FALSE!!!!

#get rid of http, https, t.co, rt, amp, single number digits, and singular letters
troll_tweets_cleaned <- troll_tweets_cleaned %>%
  filter(word != "http", 
         word != "https",
         word !=  "t.co",
         word != "rt",
         word != "amp", 
         word != "1", 
         word != "2", 
         word != "3", 
         word != "4", 
         word != "5", 
         word != "6", 
         word != "7", 
         word != "8", 
         word != "9") # you can use %in% for a list of words

5. Look at a subset of the tweets to see how often the top words appear.

troll_tweets_small <- troll_tweets_cleaned %>%
  count(word) %>%
  slice_max(order_by = n, n = 50) # 50 most occurring words

# visualize the number of times the 50 top words appear
ggplot(troll_tweets_small, 
       aes(y = fct_reorder(word,n), x = n)) +
  geom_col()

6. Sentiment Analysis

1. Get the sentiments using the “bing” parameter (which classifies words into “positive” or “negative”)
2. Report how many positive and negative words there are in the dataset. Are there more positive or negative words, and why do you think this might be?

Be sure to remove the eval=FALSE!!!!

# look at sentiment
sentiments <-get_sentiments("bing")

# assign a sentiment to each word that has one associated
troll_tweets_sentiment <- troll_tweets_cleaned %>%
  inner_join(sentiments)

# count the sentiments
troll_tweets_sentiment %>%
  count(sentiment)

There are more negative words than positive because the fake tweets are trying to tear down the competition with their tweets.

7. Using the troll_tweets_small dataset, make a wordcloud:

1. That is sized by the number of times that a word appears in the tweets
2. That is colored by sentiment (positive or negative)

Be sure to remove the eval=FALSE!!!!

# make a wordcloud where the size of the word is based on the number of times the word appears across the tweets

troll_tweets_small %>%
  with(wordcloud(word, n, max.words = 12))

# make a wordcloud colored by sentiment

troll_tweets_sentiment %>%
  count(word, sentiment, sort = TRUE) %>%
  acast(word ~ sentiment, value.var = "n", fill = 0) %>%
  comparison.cloud(colors = c("blue", "orange"),
                   max.words = 20)

Are there any words whose categorization as “positive” or “negative” surprised you?

trump as positive surprised me for sure as breaking being negative just because of the many different uses of the word breaking, like breaking news, isn’t necessarily negative.

Projects

Read the project description on the moodle page. Talk to your group members about potential topics.

Task:

Write a short paragraph about ideas you have. If you already have some data sources in mind, you can link to those, but I’m more concerned with you having a topic that you’re interested in investigating right now.

We are considering doing projects related to sports or crime, we have thought of using fifa video game data to predict biases of the game, we have also come up with the idea to cluster different baseball players into different player types, and then analyzing what team makeup makes the best team. We also have discussed creating a recommender project, similar to netflix or google searches.

“Undoing” bias

Task:

Read this tweet thread by Deb Raji who you may remember from the Coded Bias film. Write a short paragraph that discusses at least one of the misconceptions.

It was interesting to me when she brought up the point that the bias can begin anywhere in the process. All the way back to the data collection process things can be biased. When you are doing your preprocessing, the process could be biased. I think that this needs to be recognized even more and thought about while working with data.